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Exoplanets are often associated with disks of dust and debris, analogs of the Kuiper 
Belt in our solar system. These "debris disks" show a variety of non-trivial structures 
attributed to planetary perturbations and utilized to constrain the properties of the 
planets 1-3 . However, analyses of these systems have largely ignored the fact that, in- 
creasingly, debris disks are found to contain small quantities of gas 4-9 , a component 
all debris disks should contain at some level 10,11 . Several debris disks have been mea- 
sured with a dust-to-gas ratio around unity 4-9 where the effect of hydrodynamics on 
the structure of the disk cannot be ignored 12,13 . Here we report that dust-gas interac- 
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tions can produce some of the key patterns seen in debris disks that were previously 
attributed to planets. Through linear and nonlinear modeling of the hydrodynamical 
problem, we find that a robust clumping instability exists in this configuration, orga- 
nizing the dust into narrow, eccentric rings, similar to the Fomalhaut debris disk 14 . The 
hypothesis that these disks might contain planets, though thrilling, is not necessarily 
required to explain these systems. 

Disks around young stars appear to pass through an evolutionary phase when the 
disk is optically-thin and the dust-to-gas ratio £ ranges from 0.1 to 10. The nearby stars 
/3 Pictoris 5,6 ' 15-17 , HD32297 7 , 49 Ceti 4 , and HD 21997 9 , all host dust disks resembling or- 
dinary debris disks and also have stable circumstellar gas detected in molecular CO, Na 
I or other metal lines; the inferred mass of gas ranges from Lunar masses to a few Earth 
masses (see Supplementary Information, Sectl). The gas in these disks is thought to be 
produced by planetesimals or dust grains themselves, via sublimation, photodesorption 10 
or collisions 11 , processes that should occur in every debris disk at some level. 

Structures may form in these disk via a recently proposed instability 12 ' 13 . Gas drag 
causes dust in a disk to concentrate at pressure maxima 18 ; but when the disk is optically- 
thin to starlight, the gas is most likely primarily heated by the dust, via photoelectric 
heating. In this circumstance, a concentration of dust that heats the gas creates a local 
pressure maximum that in turn can cause the dust to concentrate more. The result of 
this photoelectric instability could be that the dust clumps into rings or spiral patterns or 
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other structures, that could be detected via coronographic imaging or other methods. 


Indeed, images of debris disks and transitional disks show a range of asymmetries 
and other structures that beg for explanation. Traditionally, explanations for these struc- 
tures rely on planetary perturbers - a tantalizing possibility. But so far, it has been difficult 
to prove that these patterns are clearly associated with exoplanets 19,20 . 

Previous investigations of hydrodynamical instabilities in debris disks neglected a 
crucial aspect of the dynamics: the momentum equations for the dust and gas. Equilib- 
rium terminal velocities are assumed between time-steps in the numerical solution, and 
the dust distribution is updated accordingly. The continuity equation for the gas is not 
solved, i.e., the gas distribution is assumed to be time-independent, despite heating, cool- 
ing, and drag forces. Moreover, prior investigations only considered one-dimensional 
models, which can only investigate azimuthally-symmetric, ring-like patterns. This limi- 
tation also left open the possibility that in higher dimensions, the power in the instability 
might collect in higher azimuthal wavenumbers, generating only unobservable clumps. 

We present simulations of the fully compressible problem, solving for the continu- 
ity, Navier-Stokes, and energy equations for the gas, and the momentum equation for the 
dust. Gas and dust interact dynamically via a drag force, and thermally via photoelectric 
heating. These are parametrized via a dynamical coupling time t y, and a thermal cou- 
pling time r r (see Supplementary Information, Sect 2). The simulations are performed 
with the Pencil Code 21-24 , which solves the hydrodynamics on a grid. Two numerical 
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models are presented: ( 1 ) a threedimensional box embedded in the disk that co-rotates 
with the flow at a fixed distance from the star; and ( 2 ) a twodimensional global model 
of the disk in the inertial frame. In the former the dust is treated as a fluid, with a sep- 
arate continuity equation. In the latter the dust is represented by discrete particles with 
position and velocities that are independent of the grid. 

We perform a stability analysis of the linearized system of equations, that should 
help interpret the results of the simulations (see Supplementary Information, Sect 3). We 
plot in Fig. la-c the three solutions that show linear growth, as functions of e and n = kH, 
where k is the radial wavenumber and H is the gas scale height ( H=c s / where c s 

is the sound speed, O k the Keplerian rotation frequency and 7 the adiabatic index). The 
friction time t y is assumed to be equal to \ / O k . The left and middle panels show the 
growth and damping rates. The right panels show the oscillation frequencies. There is 
no linear instability for e > 1 or n < 1. At low dust load and high wavenumber the 
three growing modes appear. The growing modes shown in Fig. la have zero oscillation 
frequency, characterizing a true instability. The two other growing solutions (Fig. lb-c) are 
overstabilities, given the associated non-zero oscillation frequencies. The pattern of larger 
growth rates at large n and low e invites to take £ = en 2 as characteristic variable, and 
explore the behavior of £ 7 > 1. The solutions in this approximation are plotted in Fig. lf-g. 
The instability (red) has growth rate ~ 0.26 Ok for all £. The overstability (yellow) reaches 
an asymptotic growth rate of Ok/ 2, at ever growing oscillation frequencies. Damped 
oscillations (blue) occur at frequency close to the epicyclic frequency. 
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Whereas the inviscid solution has growth even for very small wavelengths, viscosity 
will cap power at this regime, leading to a finite fastest growing mode (see Supplemen- 
tary Information, Sect 4), that we reproduce numerically (Fig. lh). Though there is no 
linear growth for e > 1, we show that there exists nonlinear growth for e = 1. As seen 
in Fig. li, a qualitative change in the behaviour of the system (a bifurcation) occurs when 
the noise amplitude of the initial velocity (i/ rms ) is raised enough, as expected from non- 
linear instabilities 25,26 . This is an important result to emphasize because, depending on 
the abundance of H 2 , the range of £ in debris disks spans both the linear and nonlinear 
regimes. The parameter space of r T and ly is explored in ID models in the Supplementary 
Information (Sect 5), showing robustness. 

In Fig. 2 we show the linear development and saturation of the photoelectric insta- 
bility in a vertically stratified local box of size (1 x 1 x 0.6) H and resolution 255 x 256 
x 128. The dust and gas are initialized in equilibrium (see Supplementary Information, 
sect 6). The dust-to-gas ratio is given by logs = —0.75, so that there is linear instability, 
and viscosity v = a,c s H is applied as a = 10~ 4 (where a is a dimensionless parameter 27 ). 
The initial noise is u rms / c s = 10 2 . Fig. 2a shows the dust density in the x-z plane, and 
Fig. 2b in the x-y plane, both at 100 orbits. Fig. 2c shows the ID x-dependent vertical and 
azimuthal average vs time. Via photoelectric heating, pressure maxima are generated at 
the locations where dust concentrates, that in turn attract more dust via the drag force. 
There is no hint of unstable short-wavelength (< H ) nonaxisymmetric modes: the insta- 
bility generates stripes. The simulation also shows that stratification does not quench the 
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instability. Fig. 2d shows the maximum dust density vs time, achieving saturation and 
steady state around 70 orbits. 

We consider now a 2D global model. The resulting flow, in the r — (p plane (r is ra- 
dius and (p is azimuth), is shown in Fig. 3a-c at selected snapshots. The flow develops into 
a dynamic system of narrow rings. Whereas some of the rings break into arcs, some main- 
tain axisymmetry for the whole timespan of the simulation. It is also observed that some 
arcs later reform into rings. We check that in the absense of the drag force backreaction, 
the system does not develop rings (see Supplementary Information, Sect 7). We also check 
that when the conditions for the streaming instability 24 are considered, the photoelectric 
instability dominates (see Supplementary Information, Sect 8). 

An interesting development of the model is that some of the rings start to oscil- 
late, appearing eccentric. These oscillations are epicycles in the orbital plane, with a pe- 
riod equaling the Keplerian, corresponding to the free oscillations in the right hand side 
of Fig. la-c. We check (see Supplementary Information, Sect 9) that they correspond to 
eigenvectors for which u = v, that is, gas and dust velocities coinciding. For this mode, 
the drag force and backreaction are cancelled. So, for maintaining the eccentricity, this 
mode is being selected among the other modes in the spectrum. This is naturally ex- 
pected when the dust-to-gas ratio is very high. For e 1, the gas is strongly coupled to 
the dust, canceling the gas-dust drift velocity in the same way that Tf « 1 does in the 
opposite way, by strongly coupling the dust to the gas. In this configuration, the freely 
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oscillating epicyclic modes can be selected. 


We plot in Fig. 3e one of the oscillating rings, showing that its shape is better fit by 
an ellipse (red dotted line) than a circle (black dotted line). The eccentricity is 0.03, a value 
that is close to the eccentricity found 28 for the ring around HD 61005 (e=0.045 ± 0.015). We 
also notice that some of the clumps in Fig 3 should become very bright in reflected light, 
as they have dust enhancements of an order of magnitude. In conclusion, the proposed 
photoelectric instability provides simple and plausible explanations for rings in debris 
disks, their eccentricities, and bright moving sources in reflected light. 

Recent work 29 suggests that the ring around Fomalhaut is confined by a pair of 
shepherding terrestrial mass planets, below the current detection imits. Detection of gas 
around the ring would be a way to distinguish that scenario from the one we propose. At 
present, only upper limits on the amount of gas in the Fomalhaut system exist 30 ; they are 
relatively insensitive, however, because they probed CO emission, and CO could easily 
be dissociated around this early A type star. 


7 


1. Kuchner, M.J. & Holman, M.J. The Geometry of Resonant Signatures in Debris Disks 
with Planets. Astrophys. J., 588, 1110-1120, (2003). 

2. Chiang, E., Kite, E., Kalas, P., Graham, J. & Clampin, M. Fomalhaut's Debris Disk and 
Planet: Constraining the Mass of Fomalhaut b from disk Morphology. Astrophys. /., 
693 , 734-749, (2009). 

3. Lagrange, A.-M. et al. A Giant Planet Imaged in the Disk of the Young Star /i Pictoris. 
Science, 329 , 57-60, (2010). 

4. Zuckerman, B., Forveille, T., & Kastner, J. H. Inhibition of giant-planet formation by 
rapid gas depletion around young stars. Nature, 373 , 494-496, (1995). 

5. Lagrange, A. et al. The beta Pictoris circumstellar disk. XXIV. Clues to the origin of the 
stable gas. Astron. Astrophys., 330, 1091-1108, (1998). 

6. Roberge, A., Feldman, P.D., Weinberger, A.J., Deleuil, M., & Bouret, J-C. Stabilization 
of the disk around /3 Pictoris by extremely carbon-rich gas. Nature, 441 , 724-726, (2006). 

7. Redfield, S. Gas Absorption Detected from the Edge-on Debris Disk Surrounding HD 
32297. Astrophys. J., 656, L97-L100, (2007). 

8. Maness, H. L., Fitzgerald, M. P, Paladini, R., Kalas, P, Duchene, G.,& Graham, J. R. 
CARMA Millimeter-Wave Aperture Synthesis Imaging of the HD 32297 Debris Disk. 
Astrophys. J., 686 , L25-L28, (2008). 

9. Moor, A. et al. Molecular Gas in Young Debris Disks. Astrophys. /., 740 , L7-L12, (2011). 


8 


10. Grigorieva, A., Thebault, P v Artymowicz, P., & Brandeker, A. Survival of icy grains in 
debris discs. The role of photosputtering. Astron. Astrophys., 475 , 755-764, (2007). 

11. Czechowski, A., & Mann, I. Collisional Vaporization of Dust and Production of Gas 
in the /3 Pictoris Dust Disk. Astrophys. /., 660, 1541-1555, (2007). 

12. Klahr, H. & Lin, D. N. C. Dust Distribution in Gas Disks. II. Self-induced Ring Forma- 
tion through a Clumping Instability. Astrophys. /., 632 , 1113-1121, (2005). 

13. Besla, G. & Wu Y. Formation of Narrow Dust Rings in Circumstellar Debris Disks. 
Astrophys. /., 655 , 528-540, (2007). 

14. Kalas, R, Graham, J. R., & Clampin, M. A planetary system as the origin of structure 
in Fomalhaut's dust belt. Nature, 435 , 1067-1070, (2005). 

15. Olofsson, G., Liseau, R., & Brandeker, A. Widespread Atomic Gas Emission Reveals 
the Rotation of the /3 Pictoris Disk. Astrophys. }., 563 , L77-L80, (2001). 

16. Brandeker, A., Liseau, R., Olofsson, G., & Fridlund, M. The spatial structure of the /3 
Pictoris gas disk. Astron. Astrophys., 413 , 681-691, (2004). 

17. Troutman, M.R., Hinkle, K.H., Najita, J.R., Rettig, T.W. & Brittain, S.D. Ro-vibrational 
CO Detected in the /3 Pictoris Circumstellar Disk. Astrophys. /., 738 , 12-19, (2011). 

18. Takeuchi, T. & Artymowicz, P. Dust Migration and Morphology in Optically Thin 
Circumstellar Gas Disks. Astrophys. /., 557 , 990-1006, (2001). 


9 


19. Janson, M., Carson, J. C., Lafreniere, D v Spiegel, D. S., Bent, J. R., & Wong, R Infrared 
Non-detection of Fomalhautb: Implications for the Planet Interpretation. Astrophys. /., 
747, 116-122, (2012). 

20. Currie, T., Debes, J., Rodigas, T. J., Burrows, A., Itoh, Y., Fukagawa, M., Kenyon, S. J., 
Kuchner, M., & Matsumura, S. Direct Imaging Confirmation and Characterization of 
a Dust-enshrouded Candidate Exoplanet Orbiting Fomalhaut. Astrophys. /., 760 , L32- 
L37, (2012). 

21. Brandenburg, A. & Dobler, W. Flydromagnetic turbulence in computer simulations. 
CoPhC, 147 , 471-475, (2002). 

22. Lyra, W., Johansen, A., Klahr, H., & Piskunov, N. Global magnetohydrodynamical 
models of turbulence in protoplanetary disks. I. A cylindrical potential on a Cartesian 
grid and transport of solids. Astron. Astrophys., 479 , 883-901, (2008). 

23. Lyra, W., Johansen, A., Zsom, A., Klahr, H., & Piskunov, N. Planet formation bursts 
at the borders of the dead zone in 2D numerical simulations of circumstellar disks. 
Astron. Astrophys., 497 , 869-888, (2009). 

24. Youdin, A. & Johansen, A. Protoplanetary Disk Turbulence Driven by the Stream- 
ing Instability: Linear Evolution and Numerical Methods. Astrophys. J., 662 , 613-626, 
(2007). 

25. Stuart, J. T. Nonlinear stability theory. Ann. Rev. Fluid Mech., 3 , 347-370, (1971). 


10 


26. Lesur, G. & Papaloizou, J. C. B. The subcritical baroclinic instability in local accretion 
disc models. Astron. Astrophys., 513 , 60-71, (2010). 

27. Shakura, N. I. & Sunyaev, R. A. Black holes in binary systems. Observational appear- 
ance. Astron. Astrophys., 24 , 337-355, (1973). 

28. Buenzli, E., Thalmann, C., Vigan, A., Boccaletti, A., Chauvin, G., Augereau, J. C., 
Meyer, M. R., Menard, F., Desidera, S., Messina, S., Henning, Th., Carson, J., Mon- 
tagnier, G., Beuzit, J. L., Bonavita, M., Eggenberger, A., Lagrange, A. M., Mesa, D., 
Mouillet, D., & Quanz, S. P. Dissecting the Moth: discovery of an off-centered ring in 
the HD 61005 debris disk with high-resolution imaging. Astron. Astrophys., 524 , L1-L4, 
(2010). 

29. Boley, A. C., Payne, M. J., Corder, S., Dent, W. R. F., Ford, E. B., & Shabram, M. Con- 
straining the Planetary System of Fomalhaut Using High-resolution ALMA Observa- 
tions. Astrophys. /., 750, L21-L24, (2012). 

30. Liseau, R. Molecular line observations of southern main-sequence stars with dust 
disks: alpha PS A, beta Pic, epsilon ERI and HR 4796 A. Does the low gas content of 
the beta PIC and varepsilon ERI disks hint of planets? Astron. Astrophys., 348 , 133-138, 
(1999). 

Supplementary Information is linked to the online version of the the paper at www.nature.com/ nature 

Acknowledgments The writing of this paper started at the American Museum of Natural His- 


11 


tory, with financial support by the National Science Foundation under grant no. AST10-09802, 
and completed at the Jet Propulsion Laboratory, California Institute of Technology, under a con- 
tract with the National Aeronautics and Space Administration. This research was supported by an 
allocation of advanced computing resources supported by the National Science Foundation. The 
computations were performed on the Kraken system at the National Institute for Computational 
Sciences. We acknowledge discussions with H. Latter and G. Stewart. 

Author Contributions W.L. contributed to developing the model, performed the calculations, 
and wrote the manuscript. M.K. contributed to developing the model and writing of the manuscript. 

Author information Reprints and permissions information are available at www.nature.com/ reprints. 
The authors declare no competing financial interests. Correspondence and requests should be ad- 
dressed to W.L. (wlyra@jpl.nasa.gov) or M.K. (marc.j.kuchner@nasa.gov). 


12 


Figure 1 Linear analysis of the axisymmetric modes of the photoelectric instability. 
Solutions for axisymmetric perturbations ip' = tpexp(st + ikx), where ip is a small am- 
plitude, x is the radial coordinate in the local Cartesian co-rotating frame, k is the radial 
wavenumber, t is time and s the complex frequency. Positive real s means that a per- 
turbation grows, negative s that a perturbation is damped, and imaginary s represents 
oscillations. Solutions are for a = 0, iy = 1/D K , and r r =0. a-e, The five solutions as 
functions of n = kH and e. Solutions a-c show linear growth. Growth is restricted to the 
low dust-to-gas ratio (e < 1), high wavenumber ( n > 1) region. The growing modes in b-c 
have non-zero oscillation frequencies, characterizing an overstability. Conversely, solution 
a is a true instability, d-e, Solutions that correspond to damped oscillations through most 
of the parameter space. In a small region (high dust-to-gas ratio and high frequency), 
modes are exponentially damped without oscillating, f-g, Using £ = en 2 and taking the 
limit £ » 1 allows for better visualizing the three behaviors: true instability (red), oversta- 
bility (yellow), and damped oscillations (blue). The other two solutions are the complex 
conjugate of the oscillating solutions, and not shown, h, When viscosity is considered 
(tx. = 10~ 2 in this example), power is capped at high wavenumber, leading to a finite most 
unstable wavelength. The figure shows the analytical prediction of the linear instability 
growth in this case (see Supplementary Information, Sect 4) compared to the growth rates 
measured numerically. The overall agreement is excellent. The growth rates are only very 
slightly underestimated, i, Although there is no linear instability for e = 1, growth occurs 
when the amplitude of the initial perturbation (u im s ) is increased, a hallmark of nonlinear 


13 


instability. 
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Figure 2 Growth and saturation of the photoeletric instability. In this threedimen- 
sional stratified local box with linearized Keplerian shear, the main source of heating is 
photoelectric. The equilibrium in the radial direction is between stellar gravity, Coriolis, 
and centrifugal forces. In the vertical direction the equilibrium for the gas is hydrostatic, 
between stellar gravity, pressure, and the drag force backreaction. To provide a stable 
stratification, an extra pressure p b = pci is added, where c b is a sound speed associated 
with a background temperature. For the dust, a steady state is established between grav- 
ity, diffusion, and drag force. The dust continually falls to the midplane, but is diffused 
upwards. The diffusion is applied only in z, mimicking turbulent diffusion, that is in general 
anisotropic, a, x - z cut at y = 0 at 100 orbits. The instability concentrates dust in a pre- 
ferred wavelength. The resulting structures have stable stratification, b , x-y cut at the 
midplane z = 0 at 100 orbits. No non-axisymmetric instability is observed, and the dust 
forms stripes, c, Time-evolution of the vertically and azimuthally averaged density, show- 
ing the formation of well-defined rings, d, Time evolution of the maximum dust density. 
The instability saturates at ^70 orbits in this case. The slowdown compared to the growth 
rate O k / 2 predicted in Fig. 1 is because of the use of viscosity, and the background pres- 
sure needed for the stratification. The dimensionless parameter /3 = 7 (c b /c s ) 2 measures 
the strength of this term. Panel e shows that linear instability exists as long as j 6 < 1. The 
maximum growth rates drops smoothly from Q K / 2 for /3 = 0, to zero for /3 = 1. f, The 
structure formed in the dust density at t = 50 (^8 orbits) for different values of beta. At 
moderate values of beta growth still occurs at a significant fraction of the dynamical time. 
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The run shown in panels a-d used /3=0.5. 



d Max dust density vs time 




Figure 3 Sharp eccentric rings, a-c, Snapshots of the dust density in a twodimensional 
global disk in polar coordinates, at 20, 40, and 60 orbits, respectively. The photoelectric 
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instability initially concentrates the dust axisymmetrically into rings, at a preferred wave- 
length. As the simulation proceeds, some rings maintain the axisymmetry, others break 
into arcs. Some arcs rearrange into rings at later times, such as the ones at r = 0.6 and 
r = 1.0 between snapshots b and c. Though mostly axisymmetric, some rings appear 
to oscillate, appearing off-centered or eccentric, d, We measure the azimuthal spectral 
power of the density shown in snapshot c, as a function of radius. Modes from m = 0 to 
m = 3 are shown, where m is the azimuthal wavenumber. Though the ring at r = 1.5 has 
m = 0 as the more prominent mode, we show in panel e that a circle (black dotted line) 
is not a good fit. An ellipse of eccentricity e = 0.03 (red dotted line) is a better fit, though 
still falling short of accurately describing its shape. The black and red diamonds are the 
center of the circle (the star), and the center of the ellipse (a focal distance away from the 
star). 
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Supplementary Information 


1 Gas in debris disks 

Debris disks with gas represent a regime of nebular astrophysics that has only recently been quantified 31 . The archetype 
of the class, the best studied object so far, is the disk around ft Pictoris, which we use as a reference point for our study. 


1.1 Total gas mass 

It is important to note that the total mass of gas that these debris disks have is poorly known, even for the well-studied 
disk around ft Pictoris. Debris disk gas has mainly been observed in emission lines from metal ions and CO, but the 
bulk of the gas is generally assumed to be hydrogen, a component that is difficult to measure 32 . Estimates of the 
hydrogen abundance in the /3 Pictoris disk relative to the solar value range from 1 0 3 to 1; this range translates into a 
range of total gas mass from about 8 x 1O~ 4 M0 to O.8M0. 

The dust mass is better constrained. Still, care must be taken to specify the particle size range of interest, since in 
typical grain size distributions, the larger bodies contain most of the mass, yet the smaller grains are the ones that we 
can detect. The disk should have a dust mass of O.27M0 for particles smaller than 1 cm, assuming spherical grains 
with a number distribution dn/ da, tx , where a , is the grain radius 31 . Given these numbers, £ for Pictoris would 

lie in the range of 0.3 to 300. 

These numbers roughly span the range of parameters for other debris disks with gas, with 49 Ceti a notably gas-rich 
exception. Estimates for the HD 21997 disk 9 quote a dust mass of ~ O.1M0 and a gas mass of O.35M0, corresponding 
to £ « 0.3. An upper limit to the gas mass of HD 3229 is estimated 7 at ~ O.3M0; the dust mass in the same system 8 lies 
in the range of 0.02-1 M0, yielding £ in the range «0.05-3. The 49 Ceti debris disk has about the same mass of dust as 
HD21997, but substantially more gas 4 : 13M0, corresponding to £ « 0.008. 


1.2 Mean free path 

To determine if the system may or may not be treated hydrodynamically, we estimate the mean free path 
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in actual systems. In Eq. (1) n,, is the number density of the constituent gas molecules or atoms, and <7 co ]| «2x 
10 l3 cm 2 the collisional cross section of hydrogen, taken to be representative. 

For 49 Ceti, a best fit for the CO column density of Nqo = 4 x 10 15 cm~ 2 at r =100 AU is reported 33 . The number 
density can be estimated from the column density, n g = N g / 2 H, where H is the gas scale height. For an aspect ratio 
h = H/r = 0.02, the midplane number density is Hqq = 70 cm~ 3 at 100 AU . This provides an upper limit to the 
mean free path of A m f„ < 0.4 AU. The actual mean free path is likely to be much lower, because of the presence of 
hydrogen. Assuming a hydrogen abundance up to «h 2 / n CO = 10 4 , the mean free path ranges from 4 x 10“ 5 to 0.4 AU. 
Considering typical ring structures of 10 AU, the lower limit is well within the hydrodynamical range, and even the 
upper limit implies marginal applicability. 

For ft Pictoris, observations report a column density of atomic C II of 10 16 cm~ 2 at ~ 100 AU 6 . With aspect ratio 
h= 0.2 34 , this translates into a number density of « 20 cm~ 3 and thus an upper limit on the mean free path of «2 AU. As- 
suming solar composition (»h ~ 10 3 cm~ 3 ) 34 , the mean free path should be A m fp = 0.03 AU, well in the hydrodynamic 
limit. 

For HD 21997, a simultaneous fit of the CO /= 3-2 and /=2-l lines 9 yields Hqq = 22 ± 5 cm~ 3 and «h 2 / Wco = 1000 
± 500 at 63 AU (the best fit of Hqq = 10 cm~ 3 for the /= 3-2 line assuming the canonical value Kh 2 /«co = 10 4 gives a 
poor fit for the /=2-l line). This inferred number density of «h 2 = 11 — 2.2 x 10 4 cm 3 for hydrogen places the mean 
free path at A m f„ = 0.0015 — 0.003 AU, also well within the hydrodynamic limit. 
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The situation for HD 32297 is less constrained, but limits can still be derived. The gas mass is estimated from 
Na instead of CO, with a column density of Nivj a i ~ 10 11 cm~ 2 . Assuming that the abundances are similar to those 
of /3 Pictoris 6 , and using as constrains Na measurements 16 and H I upper limits 35 , the hydrogen density should be 
NHi/Njyjai ~ 10 9 . For an aspect ratio 7z=0.1 at 100 AU 36 , this yields a number density of Hhi=3 x 10 5 cm~ 3 , and thus 
a mean free path A m f„ = 1CD 4 AU. Given the assumptions made, this is likely an overestimate. Nevertheless, given an 
integral scale of 10 AU, a revised value three orders of magnitude upward would still be in the range of applicability 
of hydrodynamics. 

1.3 Thermal time 

Another important quantity for the instability we investigate is the thermal time scale. In the model quoted 31 , the dust 
is concentrated in a ring about 100-140 AU from the star. At the peak in the dust density, the midplane gas density 
is about 10 cm -3 , the dust temperature is roughly 100 K, the gas temperature is roughly 70 K. The gas is primarily 
heated by photoelectric emission from dust grains, and primarily cooled through the C II 157.7 fim line emission and 
the total heating /cooling power is roughly 2 x 10~ 19 erg s _1 cm~ 3 . Since the specific heat of molecular hydrogen at 
70 K is roughly 1.3 x 10 8 erg g~ 4 K -1 , the thermal time scale in this model is about 0.5(nH 2 /10 cm 3 ) years. Given the 
range of possible hydrogen abundances in the disk, the range of time scales of interest corresponds to about 1 0 4 to 0.1 
orbital periods. 


2 The model 

We work primarily in the thin disk approximation, using the vertically integrated equations of hydrodynamics 
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In these equations, Eg and E d are the vertically integrated gas and dust densities, respectively; u stands for the velocity 
of the gas parcels, P is the vertically-integrated pressure, and <P is the gravitational potential. S = c v (InP — qlnUg) is 
the gas entropy, where c v is the specific heat at constant volume and 7 = c p / c v is the adiabatic index, with c p the heat 
capacity at constant pressure. T stands for the gas temperature. A tridimensional model is considered in Sect. 6. 

The dust evolves Lagrangianly according to 


d i = ~™ + fi ( 6 ) 

where x is the position of a dust particle and v its velocity. The gravitational potential is given by = — G /M , / r, 
where G is the gravitational constant, M* the stellar mass, and r the stellocentric distance. The pressure is given by 
P = DgC 2 / 7, where c s is the sound speed. 

The model is closed by specifying the drag force f d by which gas and dust interact; and T p , a simple prescription 
for the gas temperature set by photoelectric heating. These are given by 


u 

T P 


(p - «) 
T / 


(7) 

( 8 ) 


The quantities t y and t t are the dynamical and thermal coupling times between gas and dust, respectively. They 
have a radial variance to match the Keplerian rate 
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Table 1: Symbols used in this work 


Symbol 

Definition 

Description 

f 


time 

U, V 


gas and dust velocity 

p , p 


gas pressure (3D and vertically integrated) 

Pg, £p 


volume and surface gas density 

Pdr £ d 


volume and surface dust density 

S 


entropy 

Cv/ Cp 


heat capacity at constant volume and at constant pressure 

T,T p 


gas and dust temperature 

Tr,T/ 


thermal relaxation and friction time 

7 


adiabatic index 

r, (p 


radial polar coordinate, azimuth 

Cs 


sound speed 

fl 


angular frequency 

s, cr 

<7 = S / fl 

complex eigenfrequency 

k, n 

n = kH 

radial wavenumber, normalized radial wavenumber 

m 


azimuthal wavenumber 

V, Oi 

v = oic s H 

viscosity and Shakura-Sunyaev parameter 

V 

v=tx'yn 1 2 

normalized viscosity 

H 

H = c s /^yfl 

gas scale height 

h 

h = H/r 

aspect ratio 

Hd 

H d = ^DJl 

dust scale height 

£ 


dust-to-gas ratio 

Cb 

Pb = pA 

sound speed associated with background temperature 

p 

£ = 7 {c b /c s f 

dimensionless background pressure parameter 

D 


dust diffusion coefficient 

£ 

v z = -£z 

proportionality factor in dust vertical velocity 


T f = T f0 n 0 /n (9) 

t t = t to fl 0 /fl (10) 

where fl = %/ GM*/r 3 is the Keplerian angular frequency. The subscript “0" refers to the values of the quantities at an 
arbitrary reference radius Tq. The quantities Tyg and Tjq are free parameters of the model. 

Given that the thermal time is sometimes expected to be very low (10 4 orbital periods, as estimated in Sect. 1.3) 
we also run models with instantaneous thermal coupling. For these models, we skip solving the energy equation, 
and equate T = T p according to Eq. (8). The sound speed is updated accordingly. This change effectively amounts to 
choosing a new equation of state that depends on the dust density 

hm P = cv( 7-1) ToSgSa/So. (11) 

A list of the mathematical symbols used in this work, together with their definitions, is provided in Table 1. We solve 
the equations with the PENCIL CODE , which integrates the evolution equations with sixth order spatial derivatives, 
and a third order Runge-Kutta time integrator. Sixth-order hyperdissipation terms are added to Eq. (2)-Eq. (4), to 
provide extra dissipation near the grid scale 22-21 . They are needed because the high order scheme of the Pencil Code 
has little overall numerical dissipation 37 . 


3 Linear stability analysis 

We perform a linear stability analysis, that should assist on interpreting the results of the numerical simulations. To 
derive the perturbation equations, we make use of the shearing sheet and fluid approximations. The first treats the 

1 The code, including improvements done for the present work, is publicly available under a GNU open source license and can be downloaded at 

http://www.nordita.org/ software/ pencil-code 
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equations in a local, co-rotating Cartesian frame. The second greatly simplifies the treatment of solid particles by 
having a continuity equation. The 2D equations are 
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where £ = T, d /E g is the dust-to-gas ratio and 


( 12 ) 

(13) 

(14) 

(15) 

(16) 

(17) 


V w = dt + w ■ V — qQxdy (18) 

is the shear-modified advective derivative 38,39 , with q = 3/2 the Keplerian shear rate. Upon linear decomposition 
ip = ip 0 + <// and considering axis-symmetric planar wave perturbations ip' = ipexp(st + ikx), these equations become 


S ^d — 

-E do ikv x 

(19) 

SV X = 

IflVy - ^{v x - u x ) 

(20) 
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(22) 

su x = 

2f2u y - Cik(£ d + eilg) - ^ (u x - v x ) 

(23) 

SUy = 

Qu x - ~{u y — Vy) 

(24) 


where we used the instantaneous thermal coupling approximation (Eq. 11) to substitute 

c 2 

VP = + E d VE g ). (25) 

Equations (19) and (22) readily allow for reducing the system to only four equations. We substitute these in the radial 
equation for gas velocity to obtain 


crv x 


av y 
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CTUx 
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u x + £\ 
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where we also substituted the normalizations 


(26) 

(27) 

(28) 
(29) 


a = s/12 (30) 

n = kc s o/y/jfl = kH (31) 
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We now substitute iy = 1 / Q so that the dispersion relation becomes simpler yet still captures the physically interesting 
case of the most mobile dust. We solve the eigenvalue problem 


(M — Icr) • A = 0 

where A = {v x ,Vy,u x ,Uy) T , I is the unit matrix, and 

-12 10 
- 1/2 -1 0 1 

e(l — n 2 /cr) 0 —£(l + n 2 /cr) 2 

0 £ - 1/2 — £ 

The dispersion relation for this linear system is a quintic polynomial 

Ait 5 + Be 4 + Ccr 3 + Da 2 + Ecr + F = 0 

with coefficients 


A 

= 1, 

(35) 

B 

= 2e + 2, 

(36) 

C 

= £ 2 + e(n 2 + 2) + 3, 

(37) 

D 

= e 2 n 2 + £(3n 2 + 2) + 2, 

(38) 

E 

= £ 2 (2n 2 + 1) + e(3h 2 + 2) + 2, 

(39) 

F 

2 2 2 
= e n — m . 

(40) 


Eq. (34) has five solutions. Since quintics do not have a complete analytical solution, we solve it numerically at first, 
to explore the behavior of the solutions. The result is detailed in the main article, and plotted in Fig la-e. In the main 
article we also explore the limit of large n and low £. Substituting £ = en 2 as characteristic variable, and taking the limit 
£ 1, the coefficients are 


(32) 


(33) 


(34) 


Ami; B = 2; C = £; 

D = 3f; E = 3£; F = -£. (41) 

The solutions are plotted in Fig lf-g of the main article. 


4 Comparing linear theory and simulations 

From the solutions (Fig la-c of main article), we see that there is significant growth even for very small wavelengths. 
The simulations, however, will cap power at and near the grid scale. To make for a meaningful comparison, we add 
artificial Laplacian viscosity v to the gas and dust momentum equations. The extra term in Fourier space is proportional 
to vk 2 , which, using the alpha-viscosity recipe 27 v = ac s H, normalizing by Q, and substituting Eq. (31), reduces to 
v = ay/r. This enters in the coefficient matrix as diagonal terms. The new, viscous, system is therefore 

[M - I(cr + !>)]• A = 0 (42) 

We set a = 1 0 2 and solve the system numerically. Fig lh of the main article shows a comparison between the linear 
growth rates predicted by Eq. (42), and the ones we measure by solving the system in the shearing sheet with fluid 
approximation (Eqs. (12) and (17)) with the Pencil Code, applying a small amplitude perturbation (M rms = 10~ 3 c s ). 
The agreement is excellent, with the measured growth rates only very slightly systematically offset from the analytical 
prediction. 


5 Nonlinear numerical simulations 

5.1 Numerical results and Robustness 

5.1.1 The instability in one dimension 

To understand the nature of the instability, it is instructive to consider it first in one-dimensional models, since it 
allows a more thorough exploration of the parameter space. We also shut down the drag force backreaction, in order 


5 



Supplementary Figure 1: Parameter space exploration for the thermal coupling time Tjq in models with Tyo=l and 
without backreaction. a, The pressure build-up associated with too short thermal coupling time leads to modification 
of the centrifugal balance experienced by a gas parcel, b, In the case of r r =10 the pressure buildup leads to supersonic 
motion and consequent disruption of the dust fingers into a turbulent pattern, c-d. The models with larger Tjq are 
better behaved as the buildup is slow. After long enough times (not shown) they too develop shock waves that disturb 
the dust streams. The dynamics is considerably different in models including back-reaction (Suppl. Fig. 3). 



Supplementary Figure 2: Parameter space exploration for the friction time ry 0 in models with Tt- 0 =10 and without 
backreaction. There is an apparent symmetry with respect to iy = 1. a-b, As the dust couples more strongly (lower ly), 
a grain takes longer to move away from the density maxima (pressure minima) into the pressure maxima, c-d, As the 
dust decouples (larger ry) the particles take longer times to move to the local pressure maxima. The velocities never 
get too close to sonic so turbulent disruption as in the case of short thermal coupling time does not occur. 
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a ^£1=1 X,n=1 x r £l=10 x,£2=1 b 




Supplementary Figure 3: The effect of the backreaction of the drag force is illustrated above, a-b, Space-time variation 
of the dust density in models of constant friction time and varying thermal time, c-d, Models of varying friction time 
and constant thermal time. Compared to the same models without backreaction (Suppl. Fig. 1 and Suppl. Fig. 2), 
we see that the dust streams are better shepherded. This is because when radial velocities are excited, the drag force 
backreaction damps them back to zero at an e-folding time t y Eg/ Ej. The shepherding of the dust stream should occur 
as long as this quantity is smaller or similar to r T . 

to better isolate the effect of heating (but as we show later, this term will become of paramount importance in the 
phenomenology of debris disks in the presence of gas). It is illustrated in the Fig. 1 of the main article how the initial 
random overdensities quickly grow into massive accumulations of dust. We explore here the parameter space of the 
thermal coupling time t t and the dynamical coupling time Ty. The dust-to-gas ratio is set to e=l, so the regime is 
marginally nonlinear. We use Lagrangian particles for the dust. 

Models exploring the parameter space of t t J 2 are shown in Suppl. Fig. 1. The panels from left to right correspond 
to Tj, 12=1, 10, 1000, and 1000. The friction time was kept constant at Tyl2=l so that the dust is mobile. In all cases the 
dust concentrates. 

The models shown in Suppl. Fig. 1 contain some interesting features worth highlighting. As seen in the model 
of t t 12=1.0, there are instances in time, around 30 orbits, that the dust distribution rapidly migrates inwards, setting 
in another equilibrium location. In the model with r r 12=10 the dust at later times, around 50 orbits, passes from a 
thin stream to a shower, resembling the transition to turbulence seeing in, e.g., cigarette smoke. Both effects seem 
to be related to the pressure buildup. As the dust concentrates, the local temperature rises, leading to further dust 
buildup. The runaway process has to saturate at some point, and these effects seem to be manifestations of saturation. 
The pressure buildup leads to a gradual change in the centrifugal equilibrium rep 2 = rfi £ + 27“ 1 d r P. At some point the 
buildup of pressure significantly changes that relation, and the disk readjusts. In the other case, the buildup of pressure 
leads to supersonic speeds. Shock waves develop, leading to the disruption of the dust stream. These effects will be 
significantly mitigated by the inclusion of the drag force backreaction, in Sect. 5.1.3. 

5.1.2 Effect of gas drag 

The models shown in the previous section were run with a fixed friction time of Tyl2=l, representing the most mobile 
particles. We consider now models with varying Ty, keeping t t 12=10. A suite of such models is shown in Suppl. Fig. 2. 
As the dust decouples (larger t y), it takes longer for the pressure gradient to move the particles to the pressure maxima. 
As the dust couples more strongly (lower Ty), it takes longer for the dust to decouple itself from the density maxima 
(pressure minima) and into the pressure maxima. The symmetry with respect to t yl2 = 1 is striking. 

In the extreme of particle tracers ( ry = 0) there is no instability as the gas and dust cannot get separated. The heating 
leads to expansion, that carries away the dust, leading to cooling. The process is self-regulated. In the opposite extreme 
of decoupled dust Ty — > oo there is no instability either as the dust does not get pushed toward pressure maxima. 
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5.1.3 Effect of gas drag back-reaction 

We now turn to the effect of the last term in the right hand side of the momentum equation, the drag force backreaction 
from the particles onto the gas. We re-run models of r T i?=l and 10, with Tyo=l (the leftmost ones in Suppl. Fig. 1), 
but now including this term. The results are shown in the upper panels of Suppl. Fig. 3, for the dust density only. 
Comparing these models with the ones without back-reaction shown in Suppl. Fig. 1, we see that the jerks in position 
are absent, as well as the dispersion of the particles into turbulent streams as seen in the models without backreaction. 
In other words, the backreaction of the drag force has the effect of shepherding the dust streams. 

The reason is because in the other models, although the dust is forced to follow the gas, the gas is unconstrained 
by the dust. When backreaction is included, the gas follows the dust whenever the dust-to-gas ratio £ is high. When 
£ >~ 1 the backreaction partially dominating the gas motion effectively herds the dust; a linear perturbation to the 
system executes exponentially damped oscillations, as shown in Fig. 1 of the main article. This effect will become 
important, in higher dimensions, on damping non-axisymmetric modes and confining the dust concentration into 
rings instead of clumps. 

Interestingly, this effect, being proportional to Ty, should break the symmetry with respect to Ty = 1 seen in the 
models without backreaction (Suppl. Fig. 2). Indeed, this is what is seen in the lower panels of Suppl. Fig. 3, where we 
show the result of two simulations with r r =10 and backreaction included, varying t y. The left panel has r y = 0.1 and 
the right panel has Ty = 10. The simulation with shorter friction time experiences less clumping for the same thermal 
coupling time. 

The herding provided by the drag force also introduces a dependency on r T / Ty, because if t t -C Ty, the pressure 
builds much faster than the dust can respond. 


6 Three dimensions and stratification 

We assess now the impact of three-dimensionality and stratification in the instability. We need first to establish the con- 
dition of vertical hydrostatic equilibrium. We treated so far the gas-dust debris disks as a system where the only source 
of heating was photoelectric, which leads to the equation of state Eq. (11) under instantaneous thermal coupling. This 
is a suitable approximation for one-dimensional models as in the radial direction it is the centrifugal force that balances 
gravity. However, in the vertical direction, as pressure plays the central role of establishing hydrostatic equilibrium 
against stellar gravity, photoelectric heating alone cannot sustain the gas column. Extra sources of pressure must be 
considered. 


6.1 Extra pressure 

We add to the equation of state a term that embodies a background temperature T),, not associated with the dust, from, 
e.g., H 2 photodissociation, H 2 formation, and H 2 collisional de-excitation. For simplicity, we consider this background 
temperature constant. This temperature will enter the equation of state as an extra isothermal term. For reasons of 
symmetry, we can define a "background sound speed" Cj, associated with this temperature, and thus write pb = pgcfy 
The full pressure is thus 


P = Pg c b + P(Pd)- (43) 

The term p(pd) °c PgPd is the 3D equivalent to Eq. (11). We assess the effect of the extra background pressure in the 
photoelectric instability. For the axisymmetric modes, this term contributes an extra pressure force to be added to the 
u x perturbation equation. The linear system then becomes 


where 


-12 1 0 
- 1/2 -101 
e(1 — n 1 / a) 0 — (e + fi)(n 2 /cr) — e 2 

0 £ - 1/2 — £ 


£ = 7 



2 


(44) 


(45) 


is a dimensionless normalization. The dispersion relation is not much different. It is also a quintic polynomium, with 
coefficients 


A = 1 (46) 

B = 2e + 2 (47) 

C = n 2 fi + e(n 2 + 2) + e 2 + 3 (48) 

D = e(« 2 (/3 + 3) + 2) + 2n 2 /3 + £ 2 m 2 + 2 (49) 

E = £(n 2 (fi + 3) + 2) + 2h 2 /3 + e 2 (2n 2 + 1) + 2 (50) 

F = £« 2 (/3 — 1) + £ 2 n 2 . (51) 


These coefficients resume to Eqs. (34)-(40) for /3 = 0. We show in Fig 2e of the main article the maximum growth rate 
as a function of /3. The linear instability weakens as /3 increases, eventually disappearing when |S = 1. For (5 = 0.1 the 
maximum growth rate is approximately half as that for = 0, so the instability still occurs at timescales comparable to 
the dynamical time. We show the growth for different modes in Fig 2f of the main article. There are no major differences 
qualitatively. The modes for higher /3 develop similar structures at a slower pace. We conclude that linear instability 
exists as long as photoelectric heating dominates over other heating processes. 

6.2 Stratification 

We assess now the impact of stratification in the instability. We start by finding the equilibrium structure of the dust 
sub-disk in the z-direction 40,41 . 


6.3 Dust steady state 

Like in gas-rich protoplanetary disks, the dust, once decoupled, will starting settling toward the midplane of the disk, 
obeying the z-momentum equation 


V v v z = -J? 2 z (v z - u z ). 

T f 


(52) 


Since the dust is pressureless, when v z = 0 and u z = 0 equilibrium is only possible if z = 0, i.e., the dust concentrates 
in an infinitely thin midplane layer. Support against gravity comes from v z ^ 0, so the dust is not in equilibrium, but 
in a steady state, defined by setting dt = 0 in the above equation. We also assume dy = u z = 0, leading to 


v z d z v z = -Q 2 z - t y 1 v z (53) 

The ansatz v z = —£z yields the quadratic equation £ 2 — £/xy + Q 2 = 0, so 

C = (t/ 1 ± \J tj 2 ~ 4f2 2 ^ /2. (54) 

Real solutions exist for Tf < 1 /2fi. The ansatz v z = —E,z means that matter is constantly falling toward the midplane. 
The non-zero divergence would unboundedly increase the density, if not counterbalanced by some mechanism. That 
mechanism is diffusion. The continuity equation for the dust is rewritten 


Pd = ~ Pd V • v + DV 2 p d (55) 

so that a steady state is possible, in which the infall of dust towards the midplane is balanced by diffusion away from 
it, leading to the differential equation 


This equation has a general solution 


Dd 2 p d + E,zd z p d + Cp d = 0. 


(56) 


/(z) = Ciexp^-|^j 



(57) 


where erfi is the imaginary error function. Because this function diverges at infinity where the density should tend to 
zero, the solution requires C2=0. We are thus left only with a Gaussian 


with H d 


Frf = Re[/(z)] = pdo exp 



\/D/£ the dust scale height. 


(58) 
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6.4 Gas hydrostatic equilibrium 

We now derive the equilibrium condition of the pressure-supported gas column, which follows from the gas 2 -momentum 
equation 


Setting u z = 0 and substituting v z 


_ 1 dp £ , . 

V u u z = -Q 2 z (u z - v z ). 

p dz T f 

— £z, we are left with the following condition of hydrostatic equilibrium 


(59) 


1 dp 
p dz 



(60) 


Given v z = — £z, lnp^ = — z 2 /2 H 2 , and m z = 0, we calculate p s that satisfies the equilibrium condition. Substituting 
the equation of state and e = p^/ p (now a function of z), we arrive at the following equation 


(tf (61) 

where C = c 2 0 / (jpo)- The equation above can be cast in the form 

* [y'(a + be~ z2/2c ) + dze~ z2/2c ] =- z (f~* e~ z2/2c ^j , (62) 

where y = p, a = c 2 , b = Ce o c = H 2 , d = Pdo£/ T f and / = f? 2 - The various parameters are associated with different 
physical mechanisms involved in the equilibrium: a is associated with support from the background pressure, b from 
photoelectric heating, and c from diffusion; d is related to the drag force backreaction, and / to gravity. The general 
solution for Eq. (62) is 

y(z) = F(d) + K x e z2/2c (ae l2c + fe) ~ (1+c//fl) . (63) 

For no dust backreaction (d= 0) and no photoelectric heating (b = 0) it resumes to the usual isothermal stratification 

y(z) =K 2 e~f z2/2a = poe~ z2/2H s (64) 

with Ho = C(, / Q. The term containing d is 

dc 

= e(i _ x ) V’C 1 + y)~^ +x) 2 Fi(-x, (65) 

where ip = b / a g _z2 / 2c , % = fc/a = H 2 H 2 / c 2 and 2 F 1 is the hyper-geometric function. 

We show in Suppl. Fig. 4 the resulting stratification given the same dust distribution (Htf= 0.5, £=1, r^=0.5). The 
dark yellow line shows the usual isothermal stratification. The blue line shows the effect of adding the photoelectric 
heating. The reduction in gas mass is because, with the extra pressure, less mass is needed to balance the stellar gravity. 
The reduction in mass in the center stems from the higher temperature at the dust peak. The red line adds the effect 
of the drag force backreaction. It raises the gas mass profile because it behaves as extra gravity, as dust falls onto the 
midplane dragging the gas with it. We test numerically that this equilibrium is stable. 

The model shown in Fig2a-d of the main article had box length (1 x 1 x 0.6 )H, resolution 255 x 256 x 128, £ = 1 
(that is, v z = — z, which sets r y = 0.5, according to Eq. (54)). The background temperature was set as (6 = 0.5 for the 
background pressure, with c s q = 1 for the photoelectric pressure, and dust-to-gas ratio logs = —0.75, so that there is 
linear instability. The box is started with linear noise at the percent level. The dust diffusion coefficient was D = 
5 x 10~ 3 , setting the dust scale height at Hj = \Jf5JZ, « 0.07. This diffusion is applied only in the z direction. This 
anisotropic diffusion is chosen to provide support against gravity in the vertical direction, while not quenching the 
radial instability. It is also physically motivated. Because the processes that generate and/ or shape turbulence are in 
general anisotropic (stratification, rotation, magnetic fields), turbulent diffusion is in general anisotropic as well 42,43 . In 
practice, the continuity equation has a diffusion term in the right hand side, V • /, where we substitute / = DVp (that 
yields the usual isotropic diffusion) by / = (D • V ) p g , with D = (0,0, D). 

To assess the evolution of the non-axisymmetric modes, we measured the power spectrum through the 3D sim- 
ulation, which we show in Suppl. Fig. 5. Power in x, y, and z are shown in the upper, middle, and lower panel, 
respectively. Snapshots are shown at 5 (blue), 50 (red), and 100 orbits (black). The power in x shows the development 
of the instability, at k x « 60, in agreement with the spacing between the stripes seen in Fig. 2 of the main article. The 


dp 

dz 


(c 2 + Cpd) 
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Gas stratification 
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Supplementary Figure 4: Vertical stratification of gas, in the presence of different processes, for a column of dust of 
scale height H^= 0.5 (in the same units) maintained by diffusion. The parameters b and d represent the magnitudes of 
photoelectric heating and drag force backreaction, respectively. The dark yellow line is the usual isothermal Gaussian 
stratification. When photoelectric heating from the dust column is included, the center heats up and expands, leading 
to the double-peaked profiled shown as the blue line. With the extra heating, less mass is needed to balance the stellar 
gravity. The red line shows the equilibrium profile when drag backreaction is included. The dip in the center is filled 
because the backreaction causes gas to be dragged along with the settling dust, effectively behaving as extra gravity. 


a 




c 



K 


Supplementary Figure 5: Spectral power measured in the simulation shown in Fig. 2 of the main article. The power is 
shown in three snapshots, at 5, 50, and 100 orbits, a, The instability in x appears as the conspicuous spike at k x fis 60. 
Resonant instabilities at doubled wavenumbers are also apparent, b. The power in y behaves non-monotonically in 
time. It first drops from 5 to 50, then rises again from 50 to 100 back to levels close to the initial one. There is no 
evidence for nonaxisymmetric instability, as power remains concentrated in low ky throughout the simulation, c. The 
power in z remains statistically constant in shape and magnitude, which reflects the stability of the stratification. 


11 







0.00 

1.25 2.50 3.75 

Gas Density 

5.00 0.0 

2.5 5.0 7.5 

Dust Density 

10.0 


| r r .T m -TiT) '-rrn-i' q I " ri-n- rn r r n ■ i vr . . i-p-r-r-i T T i -rr| rr.T I rn 



1 ....... I ...... ...... I . ■ I ' 


- 2-1012 - 2-1012 

X X 


Supplementary Figure 6: Fifteen orbits into a version of the fiducial model shown in Fig. 3 of the main paper, but 
zvithout the backreaction of the dragforce. Whereas the fiducial model develops into rings, this model breaks into 
several clumps. The backreaction is necessary to maintain the axisymmetry of the system, by resisting the tendency of 
the gas to counter-rotate in and around the regions of high pressure created by the dust. 


difference in power at this wavenumber between 50 and 100 orbits is also in agreement with the time development of 
the instability as seen in Fig. 2c-d of the main article. 

The power in z has remained roughly constant through the simulation. The power in y (middle panel), on the 
other hand, shows a decline from the initial noise, easily seen as the power drops five orders of magnitude from 5 to 
50 orbits. After that, between 50 and 100 orbits, which coincides with the development of the instability (as seen in 
the x-power and Fig 2 of the main article), the azimuthal power is restored to the initial levels. We conclude that no 
short-wavelength non-axisymmetric instability is associated with photoelectric heating. 


7 Excluding the drag force backreaction 

We observe that the well-ordered scenario of the photoelectric instability gives way to turbulence when the drag force 
backreaction is excluded from the calculation. We show this case in Suppl. Fig. 6. This global simulation is similar to 
the fiducial case shown in the main paper in all aspects, except that the backreaction is switched off. As seen in the 
figure, axisymmetry is heavily broken, and the disk is teared apart into clumps. The drag force backreaction clearly 
has a pivotal role in maintaining axisymmetry. We encountered this before, when we noticed in Sect. 5.1.3 that the 
ID particle streams were sharper and better confined when this term was included, though the effect is by far better 
appreciated in 2D. 

This behavior stems from a disruption of geostrophy when the backreaction is introduced. When a dust overdensity 
generates a localized high pressure region, the motion initiated by the high pressure will be made to rotate by the 
influence of the Coriolis force. If backreaction is not present, the gas is free to execute this rotation, aware of the 
presence of dust only insofar as the latter influences the pressure. Nothing precludes the existence of several such 
clumps in the same orbit, and therefore azimuthal symmetry is not preferred. 

When, however, the backreaction is included, the gas cannot revolve around freely anymore. If the gas crosses the 
dust overdensity from the inside outward, it is slowed down (if the dust concentration is low) or outright barred (if the 
dust density is high). Pressed from both sides, the dust will expand sideways, i.e., in azimuth. The overall tendency is 
therefore to smear the dust azimuthally (eventually towards a predominant m = 0 mode), and in the process keeping 
the gas and dust radially segregated. 
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Supplementary Figure 7: Comparison between the streaming and photoelectric instabilities, a. The maximum dust 
density attained in the simulations, b-d, Snapshots of dust density in r-z of the models by the the end of the simulations 
at 50 orbits, e-g, The time-evolution of the ID vertically averaged dust density. We reproduce both the streaming and 
photoelectric instabilities in isolation. When both are modeled together, the photoelectric instability dominates the 
dynamics. 
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8 Comparison with streaming instability 

We have so far considered only disks for which the initial condition consist of gas in constant pressure. This choice 
was made in order to isolate the effect of photoelectric heating by shutting down the streaming instability 24,44-46 , a 
traffic -jam clumping instability that is triggered when particles migrate (or "stream") through pressure-supported sub- 
Keplerian gas. The possibility exists then, that the proposed instability is overwhelmed by the streaming instability 
when realistic density and temperature gradients are used. We now examine this case, and show that in actuality the 
opposite is true. 

We simulate the flow in the presence of a radial pressure gradient, that will trigger the streaming instability. The 
wavelength of interest is A, ? = ;jr, where // parametrizes the strength of the pressure gradient, according to 24 

_ 1 dp _ c 2 din p ^ ( c s \ 2 

^ ~ lp s Q 2 r dr 2vj. din r 

Given the definition of the aspect ratio and scale height, this can be written ijr ss h 2 r = hH. We set a global disk 
in r and z, with radii ranging from 0.4 to 2.5. Gas and dust surface densities are initialized as power laws of index 
— 1.5, which provides the streaming. The initial dust-to-gas ratio is e = 1. We choose h 2 = 0.1 and use 768 points 
in r. The characteristic wavelength is therefore excellently resolved with 32 points. We use 32 points in z, with the 
same resolution as in radius to allow for growth of the streaming instability. Viscosity is added as a = 10 -3 . Four 
simulations are performed, switching on and off the pressure gradient and photoelectric heating. The result is shown 
in Suppl. Fig. 7. 

The upper panel shows the time-series of the maximum dust density achieved in the simulation. The middle panels 
show the state of the r — z flow at 50 orbits. The lower panels show the evolution of the vertically averaged dust density 
as a function of radius and time. The control model without pressure gradient or photoelectric heating is essentially 
constant in time and not shown. 

In the simulation without heating but with the pressure gradient (lower left panels and blue line in upper panel), 
we recover the streaming instability, with clumps forming all over the grid and achieving high densities as they trap 
particles in their streaming course. When photoelectric heating is considered but not the pressure gradient (lower right 
panels and black line in upper panel), we recover the photoelectric instability, forming structures of longer wavelength 
in radius (about A = H/2 in size), and symmetric in z. When both the pressure gradient and photoelectric heating are 
included, (lower middle panels and black line in upper panel), the maximum dust density achieved is bigger (« x2) 
than in the purely photoelectric case, but yet three to four times lower than in simulation with the streaming instability 
alone. The structures formed (middle lower panels) are of similar size as those generated in the case without the 
pressure gradient (right lower panels), and also symmetric in z. Overall, the result is intuitive. When the photoelectric 
effect is at play, gas is forced to rarefy in regions of high dust concentration, disrupting the pressure gradient responsible 
for the streaming. We conclude that it is the photoelectric mechanism, not the streaming instability, that dominates the 
dynamics. 


9 Free oscillations 

9.1 Long wavelength limit 

We find in the simulations that the rings execute oscillations. The modes in question are the free, non-damped, os- 
cillations that occur through most of the parameter space of Fig. la-c of the main article. We can find these modes 
analytically by taking the long wavelength limit (n — > 0). Ignoring the n terms, the system becomes a quartic equation, 

A — \) B = 2s + 2; C = s 2 T- 2s V 3; 

D = 2s + 2; £ = s 2 + 2s + 2; F = 0; 

which can be solved exactly. The roots are 


lim a = ±z; lim u = — (s + 1) ± i (67) 

n — >0 > 0 

i.e., two solutions are free oscillations at wave frequency 12, thus constituting epicyclic oscillations. The other two 
solutions are damped oscillations at the same frequency 12, and damping time (s + 1) -1 12 -1 . The eigenvector corre- 
sponding to the cr = ±i solution represent the particular mode for which v x = u x and Vy = liy, therefore canceling the 
drag force. This epicyclic mode is low-/c and non-damped, and should show up in the simulations as free oscillations. 

We model a ID system in the instantaneous thermal coupling approximation, as used in the linear analysis, where 
the equation of state is given by Eq. (11). The time evolution of the dust and gas densities is shown in Suppl. Fig. 8. 
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Supplementary Figure 8: Model with instantaneous thermal coupling between gas and dust. The gas (a) attains much 
denser concentrations, and the dust (b) executes high frequency oscillations. At later times the dust streams in the 
outer disk break into a roughly uniform sheet of dust extending from r=1.4 to r= 2.0 


The evolution looks a hybrid between the models considered in Sect. 5. On the one hand the dust streams are very 
thin, because of the shepherding provided by the drag force backreaction. On the other hand, the infinitesimally small 
thermal time leads to rapid wave excitation, seen in the high frequency oscillation of the dust. 

We take a close look at this behavior, zooming into the dynamics of the stream. We show in Suppl. Fig. 9 the initial 
development of the instability (ss 8 orbits) in a local box ranging x = [— 2.5,2.5]H with 256 grid points, using both 
the Eulerian fluid and Lagrangian particle approaches. Viscosity is applied to the fluid and gas at a = 1 0 2 . Fluid 
and particle approaches therefore sample different Reynolds numbers (Lagrangian particles are inviscid). The upper 
panels show the linear case, that the fluid approach (upper left) can capture. The particle simulation (upper right) will 
trigger the nonlinear regime, which is seen by the high amplitude in that case. Yet a similar wavelength was excited, 
which suggests linearity. The lower panels show the marginally nonlinear regime, at e = 1. The fluid simulation (lower 
left) was quiescent in this case for linear noise, and had to be started at a noise level of 10% of the sound speed. The 
linear mode is also present, as seen in the underdense region developing in the middle of the figure, corresponding to 
the same density dip in the upper panel. This is because the dust was depleted in that region, falling below e = 1. The 
particle simulation at the nonlinear regime (lower right) has the same overall wavelength as the fluid one, which is seen 
as the streams develop in the same location. However, the Reynolds number is different (essentially inviscid), which 
as a result leads to clearer oscillatory dynamics in the particle case. These are present in the fluid case as well, albeit 
less prominently. These oscillations are responsible for the non-zero eccentricity seen in two-dimensional models. 

These modes are u = v modes, that naturally occur at high dust-to-gas ratio, as dust drags the gas along. The ques- 
tion is therefore how to excite these modes. Excitation of the high-epsilon epicyclic mode is possible in the following 
way. If the orbits of the dust in the ring should happen to become slightly eccentric, ordinarily, gas drag would damp 
the eccentricity, by slowing the grains near pericenter and pushing the grains near apocenter forward. However, when 
gas is heated by the dust, the side of an eccentric ring that is closer to the star heats the gas more, increasing the gas 
pressure on that side. The resulting gradient in the gas pressure slows the gas orbital motion on one side of the ring, 
causing the gas streamlines to become eccentric as well, and therefore canceling the drag force, sustaining the mode. 
The azimuthal pressure gradient in turn attracts more dust, creating a tendency for apsidal alignment. Notice that this 
mechanism works best when Tj <C t y, i.e., when the thermal time is short in comparison to the friction time. This way, 
the gas in eccentricity excursion has time to heat up at the periapse and establish the pressure gradient. Otherwise, the 
gas drag has enough time to damp the dust attempt at eccentricity excursion back to the circular ring. This is precisely 
what we see in Suppl. Fig. 1 and Suppl. Fig. 8, as the oscillations only appear for short thermal time. 
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Supplementary Figure 9: Time series of the ID dust density, showing the linear and marginally nonlinear development 
of the photoelectric instability in the initial t = 50 8 orbits) of the simulations, a-b, In the £ = 0.2 case there is linear 

instability, which is retrieved in both cases, c-d, In the £ = 1 case there is only nonlinear instability. The fluid simulation 
shown in the lower left panel had to be started with noise at M rms = 0.1 c s , since for linear noise it is quiescent. The 
particle simulation at £ = 1 develops epicylic oscillations (period = 1/1? = 2n), also visible in the fluid approach but 
less clear due to the higher Reynolds number. 
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